//=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|
// Classe DiscreteFourierTransform:
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
//   Para um conjunto de dados {tempo,ordenada}_n calcula a Transformada
//  Discreta de Fourier.
//   Alem disso, calcula a janela espectral.
//
//   Eh preciso informar o valor da frequencia minima e maxima onde procurar e
//  tambem o numero de frequencias.
//   Caso nao tenha a menor ideia de quais sao (voce tah numa roubada heim!),
//  assume os seguintes valores default:
//   fmin = 1.0/( Ttotal      )    --->  inverso do tempo total de observacao
//   fmax = 1.0/( 2.0*min(dt) )    --->  Pseudo frequencia de Nyquist (para
//                                      dados nao igualmente espacados) ou
//                                      a freq. de Nyquist para dados 
//                                      igualmente espacados
//   Nf = Numero de pontos dos dados
//
//
//   Os dados nao precisam estar igualmente espacados, nem ter um certo numero
//  (como potencia de dois no metodo radix2 da fft) e nem estar ordenados.
//
//   Se nao forem passadas as frequencias fmin e fmax, os dados serao ordenados
//   com o quicksorte para que elas possam ser calculadas
//
//   ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO
//     A amplitude retornada eh o coeficiente do cosseno/seno:
//
//         f(x) = A1[ cos(w1*x) + i*sin(w1*x) ] + ...
//
//     entao, o nome pode enganar. Pode-se pensar que amplitude se refere a
//    amplitude da curva, quando na verdade se refere ao coeficiente!
//
//   ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO ATENCAO
//________________________________________________________
// A.Z. - 03/07 => Criacao
//        04/07 => Retorna a fase (em radianos) da DFT tambem
//        04/07 => Pode retornar a potencia ou a amplitude da DFT
//        09/07 => Corrigi bug que fazia os vectors internos renovarem seus
//                 valores com multiplas chamadas a classe
//                 Nova definicao para amplitude.
//=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|=|
#ifndef DFT_H
#define DFT_H

#include <cmath>
#include <vector>
#include <iostream>
#include "aborta.hpp"
#include "quickSort.hpp"
#include "meusTipos.hpp"
#include "vecUtilities.hpp"
using namespace std;




class DiscreteFourierTransform{
public:
  DiscreteFourierTransform(void){};

  void push_dataF( MTvd, MTvd, MTcvd& );
  void push_dataF( MTvd, MTvd, int=-1, double=-1, double=-1 );
  void push_dataP( MTvd, MTvd, MTcvd& );
  void push_dataP( MTvd, MTvd, int=-1, double=-1, double=-1 );

  inline MTvd frequencies( void ) const {
    if( _freq.empty()    ) aborta("A DFT ainda nao foi calculada!");
    return _freq;
  }
  inline MTvd periods( void ) const {
    if( _period.empty()  ) aborta("A DFT ainda nao foi calculada!");
    return _period;
  }
  inline MTvd ReDFT(void) const{
    if( _ReDFT.empty()   ) aborta("A DFT ainda nao foi calculada!");
    return _ReDFT;
  }
  inline MTvd ImDFT(void) const{
    if( _ImDFT.empty()   ) aborta("A DFT ainda nao foi calculada!");
    return _ImDFT;
  }
  inline MTvd phase( void ) const {
    if( _phase.empty()   ) aborta("A DFT ainda nao foi calculada!");
    return _phase;
  }
  inline MTvd power( void ) const {
    if( _Power.empty()   ) aborta("A DFT ainda nao foi calculada!");
    return _Power;
  }
  inline MTvd amplitude( void ) const {
    if( _Amplitude.empty()   ) aborta("A DFT ainda nao foi calculada!");
    return _Amplitude;
  }
  inline MTvd windowPower( void ) const {
    if( _Pwindow.empty() ) aborta("A DFT ainda nao foi calculada!");
    return _Pwindow;
  }

private:
  MTvd _freq, _period, _ReDFT, _ImDFT, _phase, _Power, _Amplitude, _Pwindow;
};



//------------------------------------------------------------------------------
// DESCRICAO:
//   Calcula a transformada discreta de Fourier com os pontos igualmente
//  espacados em frequencia.
//   Pode receber o numero de frequencias e as frequencias inferior e superior.
//   Se nao receber, calcula essas coisas.
//
//   Depois manda pra outra funcao calcular a DFT.
// RECEBE:
//   t    => Abscissa (tempo)
//   y    => Ordenada
//   fmin => Frequencia minima do intervalo de frequencias da transformada
//   fmax => Frequencia maxima do intervalo de frequencias da transformada
void DiscreteFourierTransform::push_dataF( MTvd t, MTvd y, int Nf,
                                           double fmin, double fmax ){
  double daux;

  // Numero de pontos nos dados
  MTci Nd = t.size();
  if( Nd != y.size() ) aborta("Os vectors de dados nao tem o mesmo tamanho!");
  if( Nd <= 1        ) aborta("Sao necessarios pelo menos 2 pontos de dados!");

  // Ordena os dados caso precise avaliar fmin ou fmax
  if( fmin==-1 || fmax==-1 ){ 
    QuickSort qs;
    qs.referenceSort<double> (t);
    qs.sortAsTheFirst<double>(y);
  }

  // Inicia o vector com as frequencias onde buscar a peridiocidade
  if( fmin==-1 ) fmin = 1.0/(t[Nd-1]-t[0]);
  if( fmax==-1 ){
    daux=1e100;
    for( int i=1; i<Nd; i++ ){
      if( t[i]-t[i-1] < daux )
        daux = t[i]-t[i-1];
    }
    fmax = 1.0/(2.0*daux);
  }
  if( Nf<1 ) Nf = Nd;
  _freq = inivec( fmin, fmax, Nf );

  _period = MTvd(Nf);
  for( int i=0; i<Nf; i++ )
    _period[i] = 1.0/_freq[i];

  this->push_dataF(t,y,_freq);
}



//------------------------------------------------------------------------------
// DESCRICAO:
//   Calcula a transformada discreta de Fourier nos pontos de frequencia do
//  do vector freq
// RECEBE:
//   t    => Abscissa (tempo)
//   y    => Ordenada
//   freq => Vector com os pontos de frequencia
void DiscreteFourierTransform::push_dataF( MTvd t, MTvd y, MTcvd &freq ){

  MTci Nd = t.size();
  if( Nd != y.size() ) aborta("Os vectors de dados nao tem o mesmo tamanho!");
  if( Nd <= 1        ) aborta("Sao necessarios pelo menos 2 pontos de dados!");
  MTcd sqrtNd = sqrt(static_cast<double>(Nd));
  _freq = freq;
  MTci Nf = _freq.size();


  // Calcula a media das ordenadas
  double mean=y[0];
  for( int i=1; i<Nd; i++ ) mean += y[i];
  mean /= static_cast<double>(Nd);


  // Inicia os vectors
  _period    = MTvd(Nf);
  _ReDFT     = MTvd(Nf);
  _ImDFT     = MTvd(Nf);
  _Amplitude = MTvd(Nf);
  _Power     = MTvd(Nf);
  _Pwindow   = MTvd(Nf);
  _phase     = MTvd(Nf);


  double w,sum1,sum2,sum3,sum4;
  for( int k=0; k<Nf; k++ ){
    _period[k] = 1.0/_freq[k];

    w = 2*M_PI*_freq[k];
    sum1=sum2=sum3=sum4=0;

    for( int i=0; i<Nd; i++ ){
      MTcd C2  = cos( w*t[i] ), S2  = sin( w*t[i] );
      MTcd dy = y[i]-mean;
      sum1 += dy*C2;
      sum2 += dy*S2;
      sum3 += C2;
      sum4 += S2;
    }

    // Por que multiplicar por 2???? Nao entendo! Mas soh assim dah a
    // amplitude correta.
    double daux = 2.0*hypot(sum1,sum2)/static_cast<double>(Nd);
    _ReDFT[k]     = sum1/sqrtNd;
    _ImDFT[k]     = sum2/sqrtNd;
    _Amplitude[k] = daux; 
    _Power[k]     = daux*daux;
    _Pwindow[k]   = hypot(sum3,sum4)/static_cast<double>(Nd);

    // Calcula a fase e depois bota no intervalo [0,2pi), pois atan2 retorna 
    // em (-pi,pi]
    daux = atan2( sum4, sum3 );
    if( daux>=-M_PI && daux<0 ) daux += 2*M_PI;
    _phase[k] = daux;
  }
}






//------------------------------------------------------------------------------
// DESCRICAO:
//   Calcula a transformada discreta de Fourier com os pontos igualmente
//  espacados no periodo.
//   Pode receber o numero de periodos e os periodos inferior e superior.
//   Se nao receber, calcula essas coisas.
//
//   Depois manda pra outra funcao calcular a DFT.
// RECEBE:
//   t    => Abscissa (tempo)
//   y    => Ordenada
//   pmin => Periodo minimo do intervalo de periodos da transformada
//   pmax => Periodo maximo do intervalo de periodos da transformada
void DiscreteFourierTransform::push_dataP( MTvd t, MTvd y, int Np,
                                           double pmin, double pmax ){
  double daux;

  // Numero de pontos nos dados
  MTci Nd = t.size();
  if( Nd != y.size() ) aborta("Os vectors de dados nao tem o mesmo tamanho!");
  if( Nd <= 1 ) aborta("Sao necessarios pelo menos 2 pontos de dados!");

  // Ordena os dados caso precise avaliar fmin ou fmax
  if( pmin==-1 || pmax==-1 ){
    QuickSort qs;
    qs.referenceSort<double> (t);
    qs.sortAsTheFirst<double>(y);
  }


  // Inicia o vector com os periodos onde buscar a peridiocidade
  if( pmax==-1 ) pmax = t[Nd-1]-t[0];
  if( pmin==-1 ){
    daux=1e100;
    for( int i=1; i<Nd; i++ ){
      if( t[i]-t[i-1] < daux )
        daux = t[i]-t[i-1];
    }
    pmin = 2.0*daux;
  }
  if( Np<1 ) Np = Nd;
  _period = inivec( pmin, pmax, Np );

  _freq = MTvd(Np);
  for( int i=0; i<Np; i++ )
    _freq[i] = 1.0/_period[i];

  this->push_dataF(t,y,_freq);
}



//------------------------------------------------------------------------------
// DESCRICAO:
//   Calcula a transformada discreta de Fourier nos pontos de periodo do
//  do vector period
// RECEBE:
//   t      => Abscissa (tempo)
//   y      => Ordenada
//   period => Vector com os pontos de periodos
void DiscreteFourierTransform::push_dataP( MTvd t, MTvd y, MTcvd &period ){
  MTci Np=period.size();
  _period = period;
  _freq = MTvd(Np);
  for( int i=0; i<Np; i++ ) _freq[i] = 1.0/_period[i];

  this->push_dataF(t,y,_freq);
}



#endif // DFT_H
